arXiv: 1507.06484vl [math.NA] 23Jul2015 



Available online at www.sciencedirect.com 

Journal 



*1' 

ScienceDirect 


ELSEVIER 

00 ( 2015 ) f-E] 

Logo 


Numerical solution of Volterra integral equations of the first kind 

with discontinuous kernels 


Ildar Muftahov a , Aleksandr Tynda b , Denis Sidorov a,c 

a Irkutsk National Research Technical University 
b Penza State University 

c Energy Systems Institute of Russian Academy of Sciences 


Abstract 

We propose the numerical methods for solution of the weakly regular linear and nonlinear evolutionary (Volterra) integral equation 
of the first kind. The kernels of such equations have jump discontinuities along the continuous curves (endogenous delays) which 
starts at the origin. In order to linearize these equations we use the modified Newton-Kantorovich iterative process. Then for 
linear equations we propose two direct quadrature methods based on the piecewise constant and piecewise linear approximation 
of the exact solution. The accuracy of proposed numerical methods is 0(l/N) and 0(l/N 2 ) respectively. We also suggest a 
certain iterative numerical scheme enjoying the regularization properties. Furthermore, we adduce generalized numerical method 
for nonlinear equations. We employ the midpoint quadrature rule in all the cases. In conclusion we include several numerical 
examples in order to demonstrate the efficiency of proposed numerical methods. 

Keywords: Volterra integral equations, discontinuous kernels, direct quadrature method, regularization, evolving dynamical 
systems, midpoint quadrature. 


Introduction 

In this article we continue our studies of the novel class of linear Volterra (evolutionary) integral equations 
(VIE) of the first kind with piecewise continuous kernels. The solution of linear integral equations of the first kind 
is of course classical problem and has been addressed by numerious authors. But only few authors studied these 
equations in case of jump discontious kernels. In general, VIE of the first kind can be solved by reduction to equations 
of the second kind, regularization algorithms developed for Fredholm equations can be also applied as well as direct 
discretization methods. 

From the other hand, it is known that solutions of integral equations of the first kind can be unstable and this 
is a well known ill-posed problem. This is due to the fact that the Volterra operator maps the considered solution 
space into its narrow part only. Therefore, the inverse operator is not bounded. It is necessary to assess the proximity 
of the solutions and the proximity of the right-hand side using the different metrics. In addition, the proximity of the 
right-hand side should be in a stronger metric. Moreover, as shown in IJjiL solutions of the VIE can contain arbitrary 
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constants and can be unlimited as t —* 0. Here readers may also refer to 117, k! hJ j^l, where the problems of 
existence, uniqueness and asymptotic behavior of solutions of equations of this type are explored. 

Evolutionary integral equations are in the core of many mathematical models in physics, economics and 
ecology. Excellent historical overview of the results concerning the VIEs of the first kind is given by H. Brunner in 
the paper “1896 - 1996: One hundred years of Volterra integral equations of the first kind” [3y . The theory of integral 
models of evolving systems was initiated in the works of L. Kantorovich, R. Solow and V. Glushkov in the mid-20th 
Century. Here readers may refer to the papers of 12] and [12]. Such theory employs the VIEs of the first kind where 
bounds of the integration interval can be functions of time. Here readers may refer e.g. to the monograph [5). These 
models take into account the memory of a dynamical system when its past impacts its future evolution. The memory 
is implemented in the existing technological and financial structure of physical capital (equipment). The memory 
duration is determined by the age of the oldest capital unit (e.g. equipment) still employed. 

The paper fl] is devoted to the construction of iterative numerical algorithm for the systems of nonlinear 
Volterra-type equations related to the Vintage Capital Models (VCMs) |5]: 


x(t) = f H(t,r, x(t))c/t, 

y(t) 

t 

f K(t, r, x(r))dr = /(f), 
WO 


t 6 [f 0 , T), t 0 < T ^ oo, 


with unknown functions x{t) and y(t) satisfying the initial conditions: y(fo) = To < to, x(r) = <fo(f), t e (—oo, tp]. 

Numerical methods which are optimal with respect to complexity order were constructed in paper lji] for 
VIEs with certain weakly singular kernels. 

First results in studies of the Volterra equations with discontinuous kernels were formulated by G.C. Evans ||4] 
in the beginning of XX century. Results in the spectral theory of integral operators with discontinuous kernels were 
obtained by A.R Khromov in his paper [8]. Some results concerning the general approximation theory for integral 
equations with discontinuous kernels are presented in paper [|jj]. 

There are several approaches available for numerical solution of Volterra integral equations of the first kind. 
One of them is to apply classical regularizing algorithms developed for Fredholm integral equations of the first kind. 
However, the problem reduces to solving algebraic systems of equations with a full matrix, an important advantage 
of the Volterra equation is lost and there is a significant increase in arithmetic complexity of the algorithms. The 
second approach is based on a direct discretization of the initial equations. Here one may face an instability of the 
approximate solution because of errors in the initial data. The regularization properties of the direct discretization 
methods are optimal in this sense, where the discretization step is the regularization parameter associated with the 
error of the source data. However, only low-order quadrature formulas (midpoint quadrature or trapezoidal formulas) 
are suitable for approximation of the integrals. The Newton-Cotes formulas, Gregory and others (the second order 
and higher orders) generate divergent algorithms. The detailed description of regularizing direct numerical algorithms 
is described in the monograph (9(1. 
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It should be noted that it is very difficult to apply these algorithms to solve the equation ( TO in the form of 
d 1 .31 because of the kernel discontinuities (11.2b as described in Section 1. The adaptive mesh should depend on the 
curves of the jump discontinuity for each number N of divisions of the considered interval and therefore this mesh 
can not be linked to the errors in the source data. It is needed to correctly approximate the integrals. 


Below we continue our studies [13 


14] and propose two approaches for the numerical solution for Volterra 


integral equations of the first kind with piecewise continuous kernels. The first approach is a direct discretization 
based on piecewise constant and piecewise linear approximations of the exact solution (the first and the second order 
of accuracy, respectively). The second approach is based on the preliminary determination of the two acceleration 
values of the unknown function and then we employ the special regularizing iterative procedure. 

The paper is organized as follows. In Section 1, we describe the problem give some statements concern¬ 
ing the existence and uniqueness of solutions of VIEs with discontinuous kernels. Section 2 is dedicated to direct 
discretization numerical methods based on the piecewise constant and piecewise linear approximation of the exact 
solution. In Section 3, we describe the regularization method for linear first kind VIEs of this class. The modified 
Newton-Kantorovich iterative process for nonlinear VIEs is suggested in Section 4. The numerical examples are given 
in Section 5. 


1. Problem statement 

The object of our interest is the following integral equation of the first kind 

t 

f w, *»)<!*= m ,zio,n 

o 

where the kernel K(t, s ) is discontinuous along continuous curves aft), i = 1,21, and is of the form 


( 1 . 1 ) 


K(t, s) = 


K\(t, s), Q'o(f) < s < aft); 
K 2 (t, s), a\(t) < s < a 2 it); 


( 1 . 2 ) 


K n {t, s), a„-i (t) < s < a nit). 

Here ao(0 = 0, ao(0 < aq(f) < • ■ ■ < a„(f) = t, /(0) = 0. Let us assume, that the kernels Kft, s ) and the right-hand 
side f(t) in the equation ( 11.1b are continuous and sufficiently smooth functions. The functions aft) e C 1 [0. T] are not 
decrescent. Moreover 

arj(0) < a' 2 i0) < ... < a^.^O) < 1. 

Let us rewrite the equation (11.1b 


»i(/) a 2 (t) t 

K\(t, s)xis) ds + J' K 2 (t, s)xis) ds + ■ ■ ■ + J' K„it, s)xis) ds — fit), t e [0, T], 

0 a\ (r) a„_i(r) 


f 


(1.3) 
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It is to be noted that conventional Glushkov integral model of evolving systems is the special case of this equation 
where all the functions K,(t, s ) are zeros except of K n {t, s). 


2. Direct discretization 

2.1. Piecewise constant approximation 

Let us introduce the mesh nodes (not necessarily uniform) to construct the numeric solution of the equation 
(fO > on the interval [0, T] (if the unique continuous solution exists) 


0 = to < t\ < t 2 < . ■ . < tN = T, h - max(f; - f,-_i) = 0(N *). 

i=l, N 


( 2 . 1 ) 


The approximate solution is determined as the following piecewise constant function 

N 

xn( t) = Y x iSi(t), t e (0, T], 6j(t ): 


i=i 


1, 16 A,- = (ti-u til, 

0, t i A i 


( 2 . 2 ) 


with the undefined coefficients jt,-, i = 1, N. We differentiate the both parts of the equation (11.3b with respect to t to 
determine a'q = x(0) 


n , “ iW 


s) 


A(i) ds + a'At)Ki(t, aj(t))x(ai(t))- 


a;-l(0 

-a' i _ l (t)K i (t, ar I _i(f))x(cr,_i(0) )■ 


From the last expression we obtain 


xq = 


f\ 0) 


2 *i(0,0) f«;(0) - cr'_,(0)| 
1=1 


(2.3) 


Here it is assumed, that the denominator of (12.3b must be not zero. We introduce the denotation fa — f{tk), k = 
. .,N and write the initial equation in the point t - t\ to define the coefficient x\: 


a,0i) 


zf 

1=1 fa , 


Kj(t i, s)x(s) ds — fa. 


(2.4) 


tt.’-ita) 


Since at this stage the lengths of all integration intervals a,(fi) - a, \(t\) in (12.41 ) don’t exceed h then based on the 
midpoint quadrature rule we have 

/i 


xi = 


Ete(ti) - ai-MUtu soitetii!!)) 


(2.5) 


i= 1 


Let us suppose now that we have already found the values * 2 , * 3 ,..., Xk- 1 . We rewrite the equation (II.lb ; 


fro. 


s)x(s) ds = f(t) ■ 


4-1 

JV 


s)xn(s) ds 


( 2 . 6 ) 


4-i 
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and we require that the last equality hoi ds for the point t = t k 

4 4-i 

K(t k , s)x(s) ds = f k - f K(t k , s)xn(s) ds. 


Ik 

!■ 

4-i 


(2.7) 


Taking into account (12.2b we have 


4 

I 


x k K(t k , s) ds = fk 


k -1 L 

Zvf 


K(t k , s) ds. 


Thus 


x k 


k~ 1 

/t - E */ f K(t k , s) ds 

7=1 

4 

f K(t k , s) ds 
4-i 


( 2 . 8 ) 


Herewith we calculate the integrals of the form f K(t k , s) ds in (12.81) using the midpoint quadrature formulas with 

O-i 

auxiliary mesh nodes related to the curves aft) of the kernels K(t. s) for each value of N. It is easy to notice that the 
error of the method is 


£n = I WO - *iv(Ollc [0-r] = 

2.2. Piecewise linear approximation 

We suppose that the approximate solution is a piecewise linear function of the following form 


(2.9) 


N i _ x 

x N (t) = V U_, + X ‘ x ‘~ l (t - fi_i) Wo, t 6 ( 0 , T], 

\ li - li-t I 


( 2 . 10 ) 


m = 


where 

1, for t e A,- = (ti-\, ti\\ 

0, for t i A;. 

We need to determine the coefficients x,-, i = 1, N, of the approximate solution. Determining by the (12.31) the coeffi¬ 
cient xq and taking into account the equality (12.71 ) we obtain 


ik 

f 

4-i 


Xk- 1 + ——— (s - tk- 1 ) I K(t k , s) ds = 
t k - a-i 


k -1 L 

f 


Xj-i + — -— (s-tj_i))/:( 4 ,s) c/s. 

0 “ 0-i 


Thus excluding x k we have 


X k — -O-1 


'r *-l t 'r - 'r 

fk ~ x k - 1 f K(t k , s) ds — 2 X;-1 f K(t k , s) ds + f (s - tj-\)K(t k , s) ds 

4-1 7 = 1 0-1 ^ ‘ 0-1 


( 2 . 11 ) 


/ C* - t k -\)K(t k , s) ds 


4-i 
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where k — 1,2,,7V. 

We approximate the integrals in (12.1 11 1 by using the midpoint quadrature formulas based on auxiliary mesh 
nodes so that the values of the functions aftf) are a subset of the set of this mesh points at each particular value of N. 
The error of this approximation method is 

Bn = I WO - XivCOIIqo.r, = 0 j ■ (2-12) 


3. Iterative method 

Let the kernels Kft, s ) be a symmetric functions in their domains, i.e. 

Kj(t, s) = Kfs, 0, 1=1,2,..., n. 

We search the approximation solution of the equation (11.11) at the mesh (12.1b as a piecewise constant function like 
(12.2b . To do this we define initial values of the xq and x\ with the formulas (12.3b . (12.5b . We rewrite the equation (11.1b 


fro. 


s)x(s) ds - f{t) ■ 




s)xn(s) ds 


(3.1) 


and designate g(t) = f(t) - J K(t, s)xy(s) ds. To define the values Xk of the required approximation solution (12.2b we 

o 

use the following iterative process: 



f t 

> 

x (m+1 \t) = x (m \f) + y 

8(f)- f 

K(t , s)x (jn \s) ds 


K n 



m = 0 , 1 ,..., 


(3.2) 


where y is a positive regularization parameter and m is a number of the iteration. 

We define the initial approximate value x (0) (f) from the aprioristic data (if we have it) of the exact solution or 
we suppose x ({)> (t) = g(t). Obviously, if the functional sequence x < " !, (t) converge to a function x 7 (t) then that function 
satisfy (HI for all y ± 0. The values Xk, k = 2,3,... ,7V, can be defined successively as 



f h 

\ 

(m+ 1 ) (m) , 

4 = 4 + r 

g(tk)~ f 

K(tk, s)x^\s) ds 


*J 

\ h 

/ 


, hi = 0,1,.... 


(3.3) 


Herewith to calculate the integrals in the (13.3b we use midpoint quadrature or trapezoidal formulas based on auxiliary 
mesh nodes related to the curves aft) of the kernels Kit, s) for each value of N. In practice, we choose the optimal 
value of the regularization parameter y with the following condition 




n " i(?) 

zf- 


Kft, s)x ( ™\s) ds - f(t) 


min 


(3.4) 


C[o,n 


for large enough m. 
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4. Nonlinear equations 


In this section we consider the extension of CD to the case of nonlinear dependency K(t , s, x(s)): 

t 

fw,sMs))Js = m l£[0,n/(0) = 0, 

0 


(4.1) 


where 


K\(t, s)Gi(s,x(s)), f, semi, 


K(t, s, x(s)) = - 


K n (t, s)G n (s, x(s)), t,s e m n . 


(4.2) 


Here m, = {t, i|Q' i _i(f) < s ^ aft)}, ao(t) = 0, a„(t ) = t, i = 1 ,n, The functions K{, f(t), aft) have continuous 
derivatives with respect to t at t, s e mj, K,ft, 1) + 0, 0,(0) = 0, 0 < aft) < 02 (f) < • • • < a n -ft) < t. The functions 
01 (f),..., a n -ft) should increase in small neighbourhood0 < f < t at least. 

The following theorem states the existence and uniqueness conditions of solution of equation (14.11 ). The proof 
is similar with proof of the Theorem 3.2 in the monograph 11311 . 


Theorem 4.1. Let fort e [0, 7] the following conditions takes place: K,(t, .v), G,(s, x(,v)j are continuous, i — l,n, aft) 

and f(t) have continuous derivatives for t, K„(t, t) 4 0, 0 = oo(f) < oi(f) < • • • < o„_i(f) < o„(f) = t for t 6 (0, T], 

o,(0) = 0, /(0) = 0. Let the functions Gfs,x(s)) satisfy Lipschitz condition |G/(s, xi(s)) - Gfs, X 2 (s))[ < q\x i - X 2 [, 

Vxi,X 2 e R 1 , q„ + n - 1 ar^fO)|Ai, 7 (0, 0) _1 (A',(0,0) — ^,+1(0,0))|(1 + qi) < 1. Then 3r > 0 such that the equation 
1=1 

( 14.71 ) has an unique local solution in Cro T i. Furthermore, if min(f — o„_i(f)) = h > 0 then the the solution can be 

r<t<T 

constructed on [r, T] using the step method combined with successive approximations. Thereby the equation m has 
the unique global solution in C[ o,r]- 


4.1. Linearization 

In order to approximate solution of (14. Il l we introduce the nonlinear integral operator 


(F x)(t) 


n “f 

= 2 J Kft, s)G,(s, 


X(s)) ds - /(f). 


(4.3) 




The equation (14. Il l can be written in an operator form as follows: 


(F x)(f) - 0. 


(4.4) 


In order to construct an iterative numerical method to equation ( 14. 1 1 , we first linearize the operator (14.3b according to 
a modified Newton-Kantorovich scheme [6]: 


x m +i = x m - [F'(xq)] 1 (F(x m )), m — 0,1,..., 


(4.5) 
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where xo(f) is the initial approximation. Then, the approximate solution of (14.4b could be determined as the following 
limit of sequence: 

x(f) = lim x m (t). (4.6) 

m —>oo 

The derivative F'(x o) of the nonlinear operator F at the point xo is defined as follows: 

j,,, . ,. F(xo + cox) - F(xq) 

F (xo) = lim-= 

( 4 — >0 Cl) 


= lim 


„ “.(0 

imif 

' 1 aj-i(t) 


Kj(t, s) [G,(i, xo(i) + wx({)) - Gj(s, xo(i))] ds. 


Implementing the limit transition under the integral sign, we finally get: 


F’(x 0 )(t) = 


n 

-s ■ 


Ki(t, s)Gi X {s, xq(s))x(s) ds, where Gt x {s, xo(s)) = 


dG,(s, x(s)) 


dx 




Thus, we obtain the operator form of Newton-Kantorovich scheme as follows: 

F'(x 0 (t))Ax m+ i(t) = -F(x m ), Ax m +i = x m+ i - x m 


(4.7) 


(4.8) 


or in the extended form 


a t (t) 
n U 

£ f «<. 

i=i ‘4. 


alt) 

" r 

s)G ix (s, x 0 (s))Ax m+ i(s) ds = f{t) - ^ I Kit, s)Gt(s, x m {s)) ds. 

;=i ^ 


at.lt) ‘~at.lt) 

We rewrite the last equation as follows 

alt) 

’ll, s)G lx (s, x 0 (.S’))x m+ i (.v) ds = T m (t ), m = 0,1,2,... 

l at-li) 

where 


L f « 

i=l 


vjj) = m + 


n 

£ f «>. 

i= 1 


s) [G lx (s, x Q (s))x m (s) - Gj(s, X„,(i))] ds. 


(4.9) 


Equations ( 14.9b are now linear Volterra equations of the first kind with respect to the unknown function x m +\(t). Note 
that the kernels Kit, s)G, x (s, xo(.v)), i — 1 ,n, remain constant during each iteration m. Since equations ( 14.9b have the 
form (11.3b we can apply the methods suggested in Section 2 and Section 3 to solve them numerically. Thus, solving 
the equations ( 14.9b . we get a sequence of approximate functions x,„ + i(f). And then, using formula (14.6b . we obtain the 
approximate solution of (ED with an accuracy depending on m. 


4.2. The convergence theorem 

Let C[0, 7] be a Banach space of continuous functions equipped with the standard norm ||x||c[o, 7 i = max \x{t)\. 
The following theorem of convergence (based on the general theory proposed in the classical monograph [61]) for iter¬ 
ative process (14.9b takes place: 
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Theorem 4.2. Let the operator F has a continuous second derivative in the sphere Flo (||x—xo|| ^ r) and the following 
conditions hold: 

1. Equation (l4~9l > has a unique solution in [0, T]for m = 0, i.e. there exists To = [-F'(xo)] *; 

2. ||Axi|| < q\ 

3. ||T 0 F"(x)|| < L, x € fio- 

If also h - Lq < \ and 1 ~^l~ 2l ' q ^ r ^ q, then equation (14.11 ) has a unique solution x* in Qq, process ( 14.91 ) 

converges to x*, and the velocity of convergence is estimated by the inequality 

l|x*-x m || < 5(i _ Vl -2h) m+ \ m = 0,1,.... 
h 

In order to prove this theorem we show that equation (14.91 ) is uniquely solvable (including the case m = 0), i.e. 
condition 1 of the theorem holds. Then we verify the boundedness of the second derivative [F"(xo)](x)) for estimating 
the constant L in condition 3. 

One can verify that the necessary condition for the second derivative [F"(xo)](x) to be bounded is a differen¬ 
tiability of the initial approximation xq( f) as well as the functions K , with respect to second variable. 


4.3. Generalized numerical method for nonlinear equations 

In this section we offer for nonlinear weakly regular Volterra equations common numerical method based on 
using midpoint quadrature rule. 

To find numerical solution of the equation ( 14. Il l on the interval [0, T] we introduce the following mesh (the 
mesh can be non-uniform) 


0 = to < t\ < t 2 < ■ ■ ■ < tN = T, h - max(f; - f,_i) = 0(N *). 

i=l,N 

Let us search for the approximate solution of the equation |4T| as following piecewise constant function 


(4.10) 


x N (t) = ^ xfft), t e (0, T], 6ft ) = 


1, for t e A,- = (f,-_i, f,-]; 
0, for t i A j 


(4.11) 


with coefficients x,, i = 1, N are under determination. In order to find xo = x(0) we differentiate both sides of the 
eauation l4.1l with respect to t: 

n / “ iW 

fit) = J' dK 'd t ’ ^ Gfs, x(s)) ds + aft) Kft, aft))Gfaff), x(aft)))~ 

!_1 ai-\(t) 


-a\_ft)Kft, a^ft))Gfai-ft),xiai-ft)))\. 


Thereby 


f(0) = 


Ml 

1 n 


dKf 0,0) 
dt 


GfO, x(0)) ds + afO)KfO, 0)Gf0, x(0)) - a'_f0)Kf0, 0)G,(0, x(0)) . 


9 
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In the last expression the coefficient xq appears in the case of nonlinear dependency. To find the coefficient xq we use 
Van Wijngaarden-Dekker-Brent method. The implementation of this method are considered in detail in 111 ]. Let 
us make the notation fk — /(ft), k = 1,..., N. The mesh point of the mesh 14.101 which coincide with aftf) we still 


denote as v,y, i.e. aftf) e A v ... Obviously v, 7 < j for i = 0, n — 1, j = 1,1V. It is to be noted that aftf) are not always 
coincide with any mesh point. Here v, ; is used as index of the segment A v .., such as aftj) 6 A Vj . (or its right-hand 
side). Let us now assume the coefficients xo, xj,..., x*_i be known. Eauation l4. ll defined in t — 4 as 


»/(« t) 


Z f * 

i= 1 v 


(4, s)Gfs,x(s)) ds - f k . 


we can rewrite as follows: hitf) + h(tk) + ■ ■ ■ + /„(/>) = fk, where 


4(4) - 


vu-i / 

zf 


4(4, x(s)) ds + 


Q-lfc) 

/ 


K\(tk, s)Gfs, x(s)) c/s, 


, 'n-l ! 

4(4) = j " 


s)G„(s, x(s)) + z / 


1. If v p _ i x- + v P ' k , p - 2,..., n - 1, then 


*7?—l,k 

4(4) = J' K p (t k , s)G p (s, x(s)) ds + 

a p -i(tk) 


Z f 

j=Vp- U+l;. 


4,(4, s)G„(s,x(s)) 


UpOk) 


f 


K p (t k , s)G p (s, x(s)) ds + 4,(4, s)G p (s, x(s)) ds. 


2. If Vp-i,* = p = 2,..., n - 1, then 


a P (tt) 

4(4)= J' K p (t k ,s)Gp(s,x(s))ds. 

a p~l Ok) 

The number of terms in each line of the last formula depends on an array v j; , defined using the input data: functions 
aft), i = l,ra — 1, and fixed (for specific N ) mesh. Each integral term we approximate using the midpoint quadrature 
rule, e.g. 


a P 0k) 

/ / \ / a D (tk) A- t v ,_i \ 

4,(4, s)G p (s,x(s)) t/s * (a p (t k ) - 4 ,- 1 ) 4 4, -- Y~^~) Gp 


<4,(4) + 4 m _i / a p (tf) + t VpJl - 1 


-,Xat 


Moreover, on those intervals where the desired function has been already determined, we select xjy(f) (he. f ^ 4-i). 
On the rest of the intervales an unknown value x* appears in the last terms. We explicitly define it and proceed in 
the loop for k. The number of these terms is determined from the initial data v, 7 analysis. To find the coefficient xo 
we also use Van Wijngaarden-Dekker-Brent method. The maximum pointwise error of proposed numerical method 
s? = max |x(4) - x ,! (f,)I has order of G(-h). 


10 
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5. Numerical examples 


5.1. Linear equations 

Let us consider the following three problems using the uniform meshes only. We define on [0, T] the two- 
mesh difference D N for the h = jj and the h — and the order of convergence p N based on the D N as follows 


with tj = tj, i = 2 j, 


D n = max |- x h 2N (tj)l 
0<i<N,0< j<2N 


. Dn 

p N = log 2 —. 

Don 


(5.1) 


(5.2) 


We use the D N and the p N to estimate the order of convergence for plorlems with unknown exact solutions. 
Let us first address the equation 


t/3 

J< 


(1 + t + i)x(i) ds - I x{s) ds - 


t 

■/ 

1/3 


(2/+1)1 V3(2f + 3)i It 2 4 


45 


18 15 


, t g [0, 2], 


with known solution x(t ) = sj2t + 1 - 1. Tab. 15.21 1 shows computed maximum pointwise errors L v = max \x(t,) 

0<i<N 


Tab. 5.1.1 The errors for various stepsizes h of the 1st example 


Piecewise constant approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.097245 

0.037330 

0.020360 

0.013031 

0.005846 

0.003016 

0.001540 

0.000781 

Dn 

0.078452 

0.027076 

0.018686 

0.009008 

0.003580 

0.002179 

0.000852 

0.000441 

Pn 

1.534798 

0.535009 

1.052673 

1.330913 

0.716153 

1,354821 

0,949923 

- 

Piecewise linear approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

9.9445E-4 

3.7228E-4 

1.1005E-4 

2.4769E-5 

7.6136E-6 

1.7694E-6 

4.8759E-7 

1.4031E-7 

Dn 

9.3805E-4 

3.3986E-4 

1.0325E-4 

2.1979E-5 

6.7433E-6 

1.6615E-6 

4.3915E-7 

1.2764E-7 

Pn 

1.464704 

1.718816 

2.231901 

1.704639 

2.020961 

1.919694 

1.782557 

- 

Iterative method 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.114115 

0.044888 

0.022213 

0.013809 

0.006252 

0.003126 

0.001475 

0.000766 

Dn 

0.088355 

0.031906 

0.019514 

0.009391 

0.003886 

0,002224 

0,000898 

0,000448 

Pn 

1.469484 

0.709318 

1.055158 

1.272992 

0.805129 

1.308369 

1.003216 

- 


x h (ti)\, the two-mesh difference D N and the order of convergence p N in double precision arithmetic applied to problem 


1 for various values of h. 
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Let us now consider the equation 

r/8 3r/8 r 

r r r t 5 6it 4 

I (1 — t ■ s)x(s) ds + | (t + s)x(s) ds - I x(s) ds =-I- 

J J J 16384 3072 

0 r/8 3r/8 

where 1(f) = f 2 is exact solution. Tab. 15.11 2 shows r; v , /> v . /r v for various values of h. 


121f 3 
384 ’ 


f e [0, 2], 


Tab. 5.1.2 The errors for various stepsizes h of the 2nd example 


Piecewise constant approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.152263 

0.084389 

0.043314 

0.021544 

0.011043 

0.005522 

0.002759 

0.001385 

d n 

0.073562 

0.041074 

0.021769 

0.010589 

0.005520 

0.002805 

0.001423 

0.000715 

Pn 

0.840735 

0.915910 

1.039741 

0.939704 

0.976789 

0.978239 

0.992681 

- 

Piecewise linear approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

2.435 IE-3 

8.3747E-4 

1.7396E-4 

8.0309E-5 

2.7166E-5 

1.3062E-5 

6.3391E-6 

3.1053E-6 

d n 

2.4588E-3 

8.3227E-4 

1.6570E-4 

6.4038E-5 

2.9543E-5 

1.5713E-5 

8.2413E-6 

4.3645E-6 

Pn 

1.562823 

2.328444 

1.371611 

1.116090 

0.910813 

0.931088 

0.917034 

- 

Iterative method 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.152264 

0.084389 

0.043315 

0.021545 

0.011043 

0.005523 

0.002760 

0.001385 

d n 

0.073563 

0.041074 

0.021770 

0.010589 

0.005521 

0.002805 

0.001424 

0.000716 

Pn 

0.840754 

0.915884 

1.039775 

0.939564 

0.976928 

0.978051 

0.991917 

- 


Finally we demonstrate the results obtained for the equation 

r/8 r/2 3r/4 r 

(1 + t + s)x(s) ds + J "(2 + ts)x(s) ds + j' (t + .v — 1 )x(,y) c/i — 4 J' x(s) ds = 

0 r/8 1/8 3r/4 




-i- l_4 - i(16f + 69Z 2 + 15 f 3 ) - e* (r - 13f + 12) + e'(4t 2 - 16f + 28) + e%(14t + 20) - 32e 2 ' 


t e [0, 2], 

with known solution T(f) = ^--4. Tab. 15.11 3 shows , D N , p N for various h. 


5.2. Nonlinear equation 

Let us consider the following equation with known solution x(f) = t + n 
r/8 r/4 r 

J(t-s) sin x(s) ds + f t (2 cos x(s)) ds + J' (-l)(sin 2 x(s) + 1) ds = 

0 r/8 r/4 
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Tab. 5.1.3 The errors for various stepsizes h of the 3rd example 


Piecewise constant approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

s 

0.449935 

0.248336 

0.127251 

0.064714 

0.033201 

0.017015 

0.008610 

0.004350 

d n 

0.201599 

0.121084 

0.062537 

0.031512 

0.016186 

0.008404 

0.004260 

0.002150 

Pn 

0.735482 

0.953223 

0.988773 

0.961195 

0.945483 

0.980078 

0.986308 

- 

Piecewise linear approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

2.1312E-2 

8.5079E-3 

4.6961E-3 

2.0226E-3 

1.8494E-3 

7.2204E-4 

3.8049E-4 

1.5079E-4 

d n 

1.2804E-2 

5.6290E-3 

3.1706E-3 

1.7218E-3 

1.3342E-3 

5.0412E-4 

2.7926E-4 

1.2197E-4 

Pn 

1.185647 

0.828102 

0.880824 

0.367916 

1.404220 

0.852152 

1.195077 

- 

Iterative method 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.397955 

0.228031 

0.119002 

0.061311 

0.033202 

0,017016 

0,008611 

0,004350 

d n 

0.169924 

0.109029 

0.057691 

0.031513 

0.016186 

0.008405 

0.004261 

0.002151 

Pn 

0.640177 

0.918293 

0.872399 

0.961200 

0.945426 

0.980055 

0.986184 

- 


= + y cos (^) + (1 + 2 0 sin (^)-2l sin (^)- ^ sin (^) + ^ sin(2t), t 6 [0, 2], 

Tab. 15.21 1 shows computed maximum pointwise errors e n = max \x(t,) - for various h. 

0 <i<N 

Tab. 5.2.1 The errors for various stepsizes h of the 4st example. 


Piecewise constant approximation 

h 

1/32 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

1/4096 

£ 

0.596074 

0.301323 

0.130566 

0.065174 

0.029381 

0.014691 

0.007345 

0.003672 


6. Conclusion 

In this article we proposed the numerical method for solution of the novel class of weakly regular linear and 
nonlinear Volterra integral equations of the first kind. We outlined the main results for this class of equation derived 
in our previous works. The main contribution of this paper are a generic numerical methods designed for solution of 
such weakly regular equations. The direct numerical methods employe the midpoint quadrature rule and have the the 
<9(1 /N) and <9(1 /N 2 ) orders of accuracy. The illustrative examples demonstrate the efficiency of proposed methods. 
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